library(factoextra)
library(cluster)
diabetes <- read.csv("winequality-red.csv", header = TRUE, as.is = FALSE, sep = ';')
diabetes <- diabetes[, -1]
head(diabetes)
##    class glucose insulin sspg
## 1 Normal      80     356  124
## 2 Normal      97     289  117
## 3 Normal     105     319  143
## 4 Normal      90     356  199
## 5 Normal      90     323  240
## 6 Normal      86     381  157

Первичный анализ данных

Описательная статистика

summary(diabetes)
##       class       glucose       insulin            sspg      
##  Chemical:36   Min.   : 70   Min.   :  45.0   Min.   : 10.0  
##  Normal  :76   1st Qu.: 90   1st Qu.: 352.0   1st Qu.:118.0  
##  Overt   :33   Median : 97   Median : 403.0   Median :156.0  
##                Mean   :122   Mean   : 540.8   Mean   :186.1  
##                3rd Qu.:112   3rd Qu.: 558.0   3rd Qu.:221.0  
##                Max.   :353   Max.   :1568.0   Max.   :748.0

Среднее и медиана отличается у всех признаков, ожидаем несимметричные распределения.

Виды признаков

Построим matrix plot для того, чтобы увидеть особенности в наших данных.

library(lattice)
library(ggplot2)
library('GGally')
ggpairs(diabetes, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))

Прологарифмируем признаки.

diabetes_l <- transform(diabetes, glucose=log(glucose), insulin=log(insulin))
names(diabetes_l)[names(diabetes_l) == 'glucose'] <- 'log_glucose'
names(diabetes_l)[names(diabetes_l) == 'insulin'] <- 'log_insulin'
ggpairs(diabetes_l, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))

Удалим единичные outliers.

А также добавим раскраску по признаку “class”.

diabetes_lo <- diabetes_l
diabetes_lo[rownames(diabetes_lo)[diabetes_lo$log_insulin < 5],] <- NA
diabetes_lo <- na.omit(diabetes_lo)
ggpairs(diabetes_lo, title="correlogram", columns=c(2:4), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = class))

При дальнейшей кластеризации классифицирующий признак– “class” рассматриваться не будет.

library("FactoMineR")
dflo <- diabetes_lo[, -1]
res.pca <- PCA(dflo, scale.unit = TRUE, ncp = 2, graph= FALSE)
get_eigenvalue(res.pca)
##       eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.10928604        70.309535                    70.30953
## Dim.2 0.83660334        27.886778                    98.19631
## Dim.3 0.05411062         1.803687                   100.00000
fviz_pca_biplot(res.pca, habillage=diabetes_lo$class)

Первые 2 ГК описывают 98% данных.

Рассмотрим 3 метода для определения числа кластеров:

  1. «Силуэт» вычисляется на основе среднего внутрикластерного расстояния и среднего расстояния до ближайшего кластера по каждому образцу.

Для точки \(i\) из кластера \(C_{I}\):

\(a(i)={\frac{1}{|C_{I}|-1}}\sum _{j\in C_{I},i\neq j}d(i,j)\) – среднее расстояние между i-той точкой и всеми остальными точками кластера.

\(b(i)=\min _{J\neq I}{\frac {1}{|C_{J}|}}\sum _{j\in C_{J}}d(i,j)\) – наименьшее среднее расстояние между i-той точкой и всеми точками любого другого кластера.

Теперь мы определяем силуэт (значение) одной точки данных:

\(s(i) = \begin{cases} 1-a(i)/b(i), & \mbox{if } a(i) < b(i) \\ 0, & \mbox{if } a(i) = b(i) \\ b(i)/a(i)-1, & \mbox{if } a(i) > b(i) \\ \end{cases}\)

Среднее s(i) по всем точками кластера - это мера того, насколько плотно сгруппированы все точки кластера. Таким образом, среднее s(i) по всем данным всего набора данных является мерой того, насколько правильно были кластеризованы данные.

  1. Метод локтя основывается на сумме квадратов расстояний между точками и центроидами.

Рассматривается характер изменения общего внутригруппового разброса с увеличением числа групп k. Объединив все n наблюдений в одну группу, мы имеем наибольшую внутрикластерную дисперсию, которая будет снижаться до 0 при \(k \rightarrow n\). На каком-то этапе снижение этой дисперсии замедляется - на графике это происходит в точке, называемой “локтем”.

  1. Статистика разрыва сравнивает общую дисперсию внутри кластера для различных значений k с их ожидаемыми значениями для распределения без кластеризации.

Пусть \(E^∗_n{log(W^∗_k)}\) обозначает оценку средней дисперсии \(W^∗_k\), полученной бутстреп-методом, когда k кластеров образованы случайными наборами объектов из исходной выборки размером n.

Тогда статистика:

\(Gap_n(k)=E^∗_n{log(W^∗_k)}−log(W_k)\)

определяет отклонение наблюдаемой дисперсии \(W_k\) от ее ожидаемой величины при справедливости нулевой гипотезы о том, что исходные данные образуют только один кластер.

При сравнительном анализе последовательности значений \(Gap_n(k),k=2,…,K_{max}\) наибольшее значение статистики соответствует наиболее полезной группировке, дисперсия которой максимально меньше внутригрупповой дисперсии кластеров, собранных из случайных объектов исходной выборки.

fviz_nbclust(dflo, kmeans, method = "silhouette", k.max = 10) #"silhouette" (for average silhouette width), "wss" (for total within sum of square) and "gap_stat" (for gap statistics).

fviz_nbclust(dflo, kmeans, method = "wss", k.max = 10)

fviz_nbclust(dflo, kmeans, method = "gap_stat", k.max = 6)

Будем рассматривать деление на 2 и 3 класстерa.

set.seed(10)
km2 <- kmeans(dflo, centers = 2, nstart = 100)
km2$size
## [1] 117  27
km2$centers#/sapply(dflo, sd)
##   log_glucose log_insulin     sspg
## 1     4.74357    6.185327 139.9145
## 2     4.60564    6.145251 378.7037
km3 <- kmeans(dflo, centers = 3, nstart = 100) 
km3$size
## [1] 52  9 83
km3$centers#/sapply(dflo, sd)
##   log_glucose log_insulin     sspg
## 1    4.562684    6.006067 239.4615
## 2    4.640043    6.213764 534.6667
## 3    4.823253    6.281513 112.4217

Кластеры получаются неравного размера (как и при рассмотрении классифицирующего признака, но там другое соотвношение индивидов в классах).

Деление на кластеры в плоскости первых 2 ГК и на ggpairs для 2 и 3 кластеров при использовании метода k-means:

fviz_cluster(km2, data = dflo)

fviz_cluster(km3, data = dflo)

ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km2$cluster)))

ggpairs(dflo, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3$cluster)))

Кластеризируем данные с меньшим числом признаков, для этого рассмотрим первые 2 ГК (на одну меньше, чем количество кластеров).

km3_pca <- kmeans(res.pca$ind$coord, centers = 3, nstart = 100) 
km3_pca$size
## [1]  18  26 100
km3_pca$centers#/sapply(res.pca$ind$coord, sd)
##        Dim.1      Dim.2
## 1 -0.8814062  1.8586132
## 2  2.8245783  0.1323704
## 3 -0.5757372 -0.3689667
fviz_cluster(km3_pca, data = res.pca$ind$coord)

#ggpairs(data.frame(res.pca$ind$coord), title="correlogram", columns=c(1:2), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(km3_pca$cluster)))

Деление также неравномерное.

Перейдем к следующему методу кластеризации– методу иерархической кластеризации с разными правилами объединения кластеров.

library(dplyr)
d <- dist(scale(dflo), method = "euclidean")
h_average <- hclust(d, method = "average")
plot(h_average, cex = 0.7, hang = -1) #агломеративная кластеризация, групповое среднее расстояние

cut_average <- cutree(h_average, k = 3)
df_cl <- mutate(dflo, cluster = cut_average)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))

h_single <- hclust(d, method = "single")
plot(h_single, cex = 0.7, hang = -1) #агломеративная кластеризация, расстояние ближнего соседа (цепочки)

cut_single <- cutree(h_single, k = 3)
df_cl <- mutate(dflo, cluster = cut_single)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))

Тут все совсем плохо ):

res.hc <- hclust(d, method = "complete") #агломеративная кластеризация, расстояние дальнего соседа (шарики)
plot(res.hc, cex = 0.7, hang = -1)
rect.hclust(res.hc, k = 3, border = 2:5)

cut_complete <- cutree(res.hc, k = 3)
df_cl <- mutate(dflo, cluster = cut_complete)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))

Дивизионная кластеризация diana:

  1. На шаге 0 все объекты объединены в один кластер.

  2. На каждом шаге кластер делится, пока на шаге n - 1 все объекты данных не будут отделены друг от друга.

  3. На каждом шаге делится кластер, назовем его R на два кластера A и B. Первоначально A равно R и B пуст. На первом этапе мы должны переместить один объект из A в B. Для каждого объекта i из A мы вычисляем среднее несходство ко всем другим объектам A:

\(d(i,A({i}))=\frac{1}{|A|-1}\sum_{i \neq j}d(i,j)\)

Объект i’, для которого уравнение выше достигает своего максимального значения, будет перемещено, поэтому положим

\(A_{new}=A_{old} \setminus \{i'\}\)

\(B_{new}=B_{old} \bigcup\{i'\}\)

На следующих этапах мы ищем другие точки для перехода из A в B. Пока A все еще содержит более одного объекта, мы вычисляем

\(d(i,A({i}))-d(i,B)=\frac{1}{|A|-1}\sum_{j \in A, i \neq j}d(i,j) - \frac{1}{|B|}\sum_{h \in B}d(i,h)\)

для каждого объекта i из A, и мы рассматриваем объект i’‘, который максимизирует эту величину. Когда максимальное значение вышеприведенного уравнения строго положительное, мы перемещаем i’’ от A до B, а затем просматриваем \(A_{new}\). С другой стороны, когда максимальное значение разности отрицательное или 0, мы останавливаем процесс, и разделение R на A и B завершено.

На каждом шаге делящего алгоритма мы также должны решить, какой кластер нужно разбить. Для этого вычислим диаметр

\(diam(Q)=max_{j \in Q, h \in Q} d(j,h)\)

для каждого кластера Q, доступного после предыдущего шага, и выбираем кластер, с наибольшим диаметром.

divisive.clust <- diana(as.matrix(d), diss = TRUE, keep.diss = TRUE) #дивизионная кластеризация
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)

cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))

Многие методы показывают похожие результаты. Если сравнивать с исходной классификацией, то методы кластеризации разбивают данные на другие группы.

divisive.clust <- diana(as.matrix(dist(scale(res.pca$ind$coord))), diss = TRUE, keep.diss = TRUE)
plot(divisive.clust, main = "Divisive", which.plots = 2, hang = -1)

cut_div <- cutree(divisive.clust, k = 3)
df_cl <- mutate(dflo, cluster = cut_div)
ggpairs(df_cl, title="correlogram", columns=c(1:3), upper = list(continuous = "points"), diag = list(continuous = "barDiag"), mapping=ggplot2::aes(colour = factor(cluster)))

Не слишком хорошо сработал метод понижения числа признаков.

library(mclust)
mc <- mclust::Mclust(dflo)
summary(mc) #первый символ относится к объему, второй - к форме, а третий - к ориентации. (E-- равный; V-- переменный; I-- расположение относительно осей координат)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -701.7836 144 29 -1547.692 -1559.786
## 
## Clustering table:
##  1  2  3 
## 33 83 28
head(mc$z) #каждому наблюдению назначается кластер с максимальной оцененной вероятностью.
##           [,1]      [,2]         [,3]
## 1 4.352532e-03 0.9956475 1.496784e-14
## 2 1.344389e-07 0.9999999 8.739538e-39
## 3 8.770285e-07 0.9999991 4.940825e-32
## 4 2.443657e-03 0.9975563 2.713471e-15
## 5 5.391036e-04 0.9994609 2.013616e-19
## 6 9.513377e-03 0.9904866 7.771781e-12
par(mfrow = c(1, 2))
plot(mc, "classification")

plot(mc, "uncertainty")

fviz_cluster(mc, data = dflo)

Результаты похожи на исходную классификацию.

Максимальный BIC (среди всех моделей) у модели VVV (но у неё же и максимальное число параметров).

Попробуем выбрать другую модель:

plot(mc, "BIC") #Зависимость критерия BIC от числа кластеров для различных вариантов параметризации ковариационной матрицы

apply(mc$BIC, 2, which.max)
## EII VII EEI VEI EVI VVI EEE VEE EVE VVE EEV VEV EVV VVV 
##   8   9   6   5   4   4   7   6   4   4   3   3   3   3

Кроме первых двух моделей BIC у остальных моделей близок, возьмем модель, кластеризующую данные на 3 кластера, но имеющюю меньше параметров.

mclustModelNames('EVV')
## $model
## [1] "EVV"
## 
## $type
## [1] "ellipsoidal, equal volume"
mc1 <- mclust::Mclust(dflo, modelNames = 'EVV')
summary(mc1)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust EVV (ellipsoidal, equal volume) model with 3 components: 
## 
##  log-likelihood   n df       BIC       ICL
##       -707.3937 144 27 -1548.972 -1560.778
## 
## Clustering table:
##  1  2  3 
## 25 88 31
head(mc1$z) 
##           [,1]      [,2]         [,3]
## 1 1.830318e-04 0.9998170 1.069044e-15
## 2 9.980166e-16 1.0000000 1.527052e-41
## 3 9.212325e-14 1.0000000 1.761841e-34
## 4 6.984441e-05 0.9999302 2.128655e-16
## 5 2.713809e-07 0.9999997 7.528674e-21
## 6 1.465955e-03 0.9985340 8.978124e-13
par(mfrow = c(1, 2))
plot(mc1, "classification")

plot(mc1, "uncertainty")

fviz_cluster(mc1, data = dflo)

Результаты немного отличаются, но все еще похожи на исходную классификацию.